{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 7: Interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Robert Johansson\n",
    "\n",
    "Source code listings for [Numerical Python - Scientific Computing and Data Science Applications with Numpy, SciPy and Matplotlib](https://www.apress.com/us/book/9781484242452) (ISBN 978-1-484242-45-2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib as mpl\n",
    "mpl.rcParams[\"font.family\"] = \"serif\"\n",
    "mpl.rcParams[\"font.size\"] = \"12\"\n",
    "\n",
    "mpl.rcParams['mathtext.fontset'] = 'stix'\n",
    "mpl.rcParams['font.family'] = 'serif'\n",
    "mpl.rcParams['font.sans-serif'] = 'stix'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from numpy import polynomial as P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from scipy import interpolate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from scipy import linalg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Polynomials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "p1 = P.Polynomial([1,2,3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "p1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "p2 = P.Polynomial.fromroots([-1, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "p2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "p1.roots()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "p2.roots()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "p1.coef"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "p1.domain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "p1.window"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p1(np.array([1.5, 2.5, 3.5]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# p1([1.5, 2.5, 3.5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "p1+p2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "p2 / 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "p1 = P.Polynomial.fromroots([1, 2, 3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "p1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "p2 = P.Polynomial.fromroots([2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "p2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "p3 = p1 // p2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "p3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "p3.roots()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "p2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "c1 = P.Chebyshev([1, 2, 3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "c1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "c1.roots()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "c = P.Chebyshev.fromroots([-1, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "l = P.Legendre.fromroots([-1, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "l"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "c(np.array([0.5, 1.5, 2.5]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "l(np.array([0.5, 1.5, 2.5]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Polynomial interpolation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "x = np.array([1, 2, 3, 4])\n",
    "y = np.array([1, 3, 5, 4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "deg = len(x) - 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "A = P.polynomial.polyvander(x, deg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "c = linalg.solve(A, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f1 = P.Polynomial(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f1(2.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "A = P.chebyshev.chebvander(x, deg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "c = linalg.solve(A, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f2 = P.Chebyshev(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f2(2.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "xx = np.linspace(x.min(), x.max(), 100)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(8, 4))\n",
    "\n",
    "ax.plot(xx, f1(xx), 'b', lw=2, label='Power basis interp.')\n",
    "ax.plot(xx, f2(xx), 'r--', lw=2, label='Chebyshev basis interp.')\n",
    "ax.scatter(x, y, label='data points')\n",
    "\n",
    "ax.legend(loc=4)\n",
    "ax.set_xticks(x)\n",
    "ax.set_ylabel(r\"$y$\", fontsize=18)\n",
    "ax.set_xlabel(r\"$x$\", fontsize=18)\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig('ch7-polynomial-interpolation.pdf');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f1b = P.Polynomial.fit(x, y, deg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f1b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f2b = P.Chebyshev.fit(x, y, deg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f2b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "np.linalg.cond(P.chebyshev.chebvander(x, deg))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "np.linalg.cond(P.chebyshev.chebvander((2*x-5)/3.0, deg))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "(2 * x - 5)/3.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f1 = P.Polynomial.fit(x, y, 1)\n",
    "f2 = P.Polynomial.fit(x, y, 2)\n",
    "f3 = P.Polynomial.fit(x, y, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "xx = np.linspace(x.min(), x.max(), 100)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(8, 4))\n",
    "\n",
    "ax.plot(xx, f1(xx), 'r', lw=2, label='1st order')\n",
    "ax.plot(xx, f2(xx), 'g', lw=2, label='2nd order')\n",
    "ax.plot(xx, f3(xx), 'b', lw=2, label='3rd order')\n",
    "ax.scatter(x, y, label='data points')\n",
    "\n",
    "ax.legend(loc=4)\n",
    "ax.set_xticks(x)\n",
    "ax.set_ylabel(r\"$y$\", fontsize=18)\n",
    "ax.set_xlabel(r\"$x$\", fontsize=18);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Runge problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "def runge(x):\n",
    "    return 1/(1 + 25 * x**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "def runge_interpolate(n):\n",
    "    x = np.linspace(-1, 1, n + 1)\n",
    "    p = P.Polynomial.fit(x, runge(x), deg=n)\n",
    "    return x, p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "xx = np.linspace(-1, 1, 250)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(8, 4))\n",
    "\n",
    "ax.plot(xx, runge(xx), 'k', lw=2, label=\"Runge's function\")\n",
    "\n",
    "n = 13\n",
    "x, p = runge_interpolate(n)\n",
    "ax.plot(x, runge(x), 'ro')\n",
    "ax.plot(xx, p(xx), 'r', label='interp. order %d' % n)\n",
    "\n",
    "n = 14\n",
    "x, p = runge_interpolate(n)\n",
    "ax.plot(x, runge(x), 'go')\n",
    "ax.plot(xx, p(xx), 'g', label='interp. order %d' % n)\n",
    "\n",
    "ax.legend(loc=8)\n",
    "ax.set_xlim(-1.1, 1.1)\n",
    "ax.set_ylim(-1, 2)\n",
    "ax.set_xticks([-1, -0.5, 0, 0.5, 1])\n",
    "ax.set_ylabel(r\"$y$\", fontsize=18)\n",
    "ax.set_xlabel(r\"$x$\", fontsize=18)\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig('ch7-polynomial-interpolation-runge.pdf');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spline interpolation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "x = np.linspace(-1, 1, 11)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "y = runge(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f = interpolate.interp1d(x, y, kind=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "xx = np.linspace(-1, 1, 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 4))\n",
    "\n",
    "ax.plot(xx, runge(xx), 'k', lw=1, label=\"Runge's function\")\n",
    "ax.plot(x, y, 'ro', label='sample points')\n",
    "ax.plot(xx, f(xx), 'r--', lw=2, label='spline order 3')\n",
    "\n",
    "ax.legend()\n",
    "ax.set_ylim(0, 1.1)\n",
    "ax.set_xticks([-1, -0.5, 0, 0.5, 1])\n",
    "ax.set_ylabel(r\"$y$\", fontsize=18)\n",
    "ax.set_xlabel(r\"$x$\", fontsize=18)\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig('ch7-spline-interpolation-runge.pdf');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "x = np.array([0, 1, 2, 3, 4, 5, 6, 7])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "y = np.array([3, 4, 3.5, 2, 1, 1.5, 1.25, 0.9])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "xx = np.linspace(x.min(), x.max(), 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 4))\n",
    "\n",
    "ax.scatter(x, y)\n",
    "\n",
    "for n in [1, 2, 3, 5]:\n",
    "    f = interpolate.interp1d(x, y, kind=n)\n",
    "    ax.plot(xx, f(xx), label='order %d' % n)\n",
    "\n",
    "ax.legend()\n",
    "ax.set_ylabel(r\"$y$\", fontsize=18)\n",
    "ax.set_xlabel(r\"$x$\", fontsize=18)\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig('ch7-spline-interpolation-orders.pdf');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multivariate interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Regular grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "x = y = np.linspace(-2, 2, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "def f(x, y):\n",
    "    return np.exp(-(x + .5)**2 - 2*(y + .5)**2) - np.exp(-(x - .5)**2 - 2*(y - .5)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "X, Y = np.meshgrid(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# simulate noisy data at fixed grid points X, Y\n",
    "Z = f(X, Y) + 0.05 * np.random.randn(*X.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f_interp = interpolate.interp2d(x, y, Z, kind='cubic')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "xx = yy = np.linspace(x.min(), x.max(), 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "ZZi = f_interp(xx, yy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "XX, YY = np.meshgrid(xx, yy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(12, 5))\n",
    "\n",
    "c = axes[0].contourf(XX, YY, f(XX, YY), 15, cmap=plt.cm.RdBu)\n",
    "axes[0].set_xlabel(r\"$x$\", fontsize=20)\n",
    "axes[0].set_ylabel(r\"$y$\", fontsize=20)\n",
    "axes[0].set_title(\"exact / high sampling\")\n",
    "cb = fig.colorbar(c, ax=axes[0])\n",
    "cb.set_label(r\"$z$\", fontsize=20)\n",
    "\n",
    "c = axes[1].contourf(XX, YY, ZZi, 15, cmap=plt.cm.RdBu)\n",
    "axes[1].set_ylim(-2.1, 2.1)\n",
    "axes[1].set_xlim(-2.1, 2.1)\n",
    "axes[1].set_xlabel(r\"$x$\", fontsize=20)\n",
    "axes[1].set_ylabel(r\"$y$\", fontsize=20)\n",
    "axes[1].scatter(X, Y, marker='x', color='k')\n",
    "axes[1].set_title(\"interpolation of noisy data / low sampling\")\n",
    "cb = fig.colorbar(c, ax=axes[1])\n",
    "cb.set_label(r\"$z$\", fontsize=20)\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig('ch7-multivariate-interpolation-regular-grid.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(6, 5))\n",
    "\n",
    "c = ax.contourf(XX, YY, ZZi, 15, cmap=plt.cm.RdBu)\n",
    "ax.set_ylim(-2.1, 2.1)\n",
    "ax.set_xlim(-2.1, 2.1)\n",
    "ax.set_xlabel(r\"$x$\", fontsize=20)\n",
    "ax.set_ylabel(r\"$y$\", fontsize=20)\n",
    "ax.scatter(X, Y, marker='x', color='k')\n",
    "cb = fig.colorbar(c, ax=ax)\n",
    "cb.set_label(r\"$z$\", fontsize=20)\n",
    "\n",
    "fig.tight_layout()\n",
    "#fig.savefig('ch7-multivariate-interpolation-regular-grid.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Irregular grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "np.random.seed(115925231)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "x = y = np.linspace(-1, 1, 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "X, Y = np.meshgrid(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "def f(x, y):\n",
    "    return np.exp(-x**2 - y**2) * np.cos(4*x) * np.sin(6*y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "Z = f(X, Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "N = 500"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "xdata = np.random.uniform(-1, 1, N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "ydata = np.random.uniform(-1, 1, N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "zdata = f(xdata, ydata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "c = ax.contourf(X, Y, Z, 15, cmap=plt.cm.RdBu);\n",
    "ax.scatter(xdata, ydata, marker='.', color='black')\n",
    "ax.set_ylim(-1,1)\n",
    "ax.set_xlim(-1,1)\n",
    "ax.set_xlabel(r\"$x$\", fontsize=20)\n",
    "ax.set_ylabel(r\"$y$\", fontsize=20)\n",
    "\n",
    "cb = fig.colorbar(c, ax=ax)\n",
    "cb.set_label(r\"$z$\", fontsize=20)\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig('ch7-multivariate-interpolation-exact.pdf');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "def z_interpolate(xdata, ydata, zdata):\n",
    "    Zi_0 = interpolate.griddata((xdata, ydata), zdata, (X, Y), method='nearest')\n",
    "    Zi_1 = interpolate.griddata((xdata, ydata), zdata, (X, Y), method='linear')\n",
    "    Zi_3 = interpolate.griddata((xdata, ydata), zdata, (X, Y), method='cubic')\n",
    "    return Zi_0, Zi_1, Zi_3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(3, 3, figsize=(12, 12), sharex=True, sharey=True)\n",
    "\n",
    "n_vec = [50, 150, 500]\n",
    "\n",
    "for idx, n in enumerate(n_vec):\n",
    "    Zi_0, Zi_1, Zi_3 = z_interpolate(xdata[:n], ydata[:n], zdata[:n])\n",
    "    axes[idx, 0].contourf(X, Y, Zi_0, 15, cmap=plt.cm.RdBu)\n",
    "    axes[idx, 0].set_ylabel(\"%d data points\\ny\" % n, fontsize=16)\n",
    "    axes[idx, 0].set_title(\"nearest\", fontsize=16)\n",
    "    axes[idx, 1].contourf(X, Y, Zi_1, 15, cmap=plt.cm.RdBu)\n",
    "    axes[idx, 1].set_title(\"linear\", fontsize=16)\n",
    "    axes[idx, 2].contourf(X, Y, Zi_3, 15, cmap=plt.cm.RdBu)\n",
    "    axes[idx, 2].set_title(\"cubic\", fontsize=16)\n",
    "\n",
    "for m in range(len(n_vec)):\n",
    "    axes[2, m].set_xlabel(\"x\", fontsize=16)\n",
    "    \n",
    "fig.tight_layout()\n",
    "fig.savefig('ch7-multivariate-interpolation-interp.pdf');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Versions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%reload_ext version_information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%version_information scipy, numpy, matplotlib"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "npbook_py310",
   "language": "python",
   "name": "npbook_py310"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
